Effectively computing transition patterns with privacy-preserved trajectory datasets

Recent advances in positioning techniques, along with the widespread use of mobile devices, make it easier to monitor and collect user trajectory information during their daily activities. An ever-growing abundance of data about trajectories of individual users paves the way for various applications that utilize user mobility information. One of the most common analysis tasks in these new applications is to extract the sequential transition patterns between two consecutive timestamps from a collection of trajectories. Such patterns have been widely exploited in diverse applications to predict and recommend next user locations based on the current position. Thus, in this paper, we explore the computation of the transition patterns, especially with a trajectory dataset collected using differential privacy, which is a de facto standard for privacy-preserving data collection and processing. Specifically, the proposed scheme relies on geo-indistinguishability, which is a variant of the well-known differential privacy, to collect trajectory data from users in a privacy-preserving manner, and exploits the functionality of the expectation-maximization algorithm to precisely estimate hidden transition patterns based on perturbed trajectory datasets collected under geo-indistinguishability. Experimental results using real trajectory datasets confirm that a good estimation of transition pattern can be achieved with the proposed method.


Introduction
Recent advances in indoor and outdoor positioning techniques make it easier to monitor and collect user trajectory information during daily activities. An ever-growing abundance of data about trajectories of individual users paves the way for various applications [1][2][3][4]. However, similar to any other novel applications, these applications that leverage large-scale user trajectory data suffer from new problems as well. This is because the personal mobility information of individuals usually contains sensitive information, such as home and workplace addresses, hospital visit records, and political affiliation, which they do not want to disclose [5,6]. Thus, most users are not comfortable to provide their trajectory information to those applications.
In the past decade, extensive studies have been conducted to protect user privacy when collecting and publishing trajectories, including spatial cloaking [7], mix-zone approach [ and encryption-based scheme [9]. Recently, as differential privacy has emerged as a de facto standard for privacy-preserving data processing, there have been several efforts to apply differential privacy to the collection and publishing of trajectory data [10][11][12][13]. These privacy-preserving methods can alleviate, to some extent, users' concerns about privacy exposure when providing their trajectory data to external service providers. However, from the service provider viewpoint, such privacy-preserving schemes lead to severe degradation in the utility of the collected dataset, eventually resulting in a loss in the quality of service. One of the most common analysis tasks in applications that utilize user mobility information is to extract the sequential transition patterns between two consecutive timestamps, which corresponds to the probability that a user located at one location in the current timestamp moves to another location in the next timestamp. Such transition patterns have been widely exploited by diverse applications when predicting and recommending the next user location based on their current position. For example, the notion of first-order Markov chain, which models the sequential transition pattern between two consecutive point-ofinterests (POIs), is used to recommend the next POI candidates to users in [2]. In [1], the sequential transition influences are integrated into the matrix factorization algorithm to enhance the accuracy of next location recommendations. In [14], transition probability matrices constructing from sequential rules are used to make location prediction for location-based-service users.
Computing transition probabilities from a collection of true trajectories is not difficult. However, the same cannot be said for the dataset consisting of privacy-preserved trajectory data. Let us consider a motivational example in Fig 1. As shown in the figure, users agree to contribute their trajectories to service providers for analysis purposes. However, owing to privacy concerns, they provide perturbed trajectories, which are obtained using geo-indistinguishability (GeoInd) [10] that is a variant of differential privacy, instead of true trajectories. During this process, it is guaranteed that individual users' true location trajectories are not disclosed to the outside of the users' devices, thus, protecting the location privacy of users. However, such a privacy-preserving mechanism leads to a loss in the utility of collected dataset. Subsequently, the resulting sequential transition patterns extracted from such dataset may have a lower utility.
To address this problem, in this paper, we explore the computation of the sequential transition patterns (i.e., transition probability between two locations) with a perturbed trajectory dataset collected using GeoInd. There have been several efforts to compute aggregate statistics from users' trajectory datasets in a privacy-preserving manner using GeoInd. Most of these works are dedicated to the estimation of population density distribution, which corresponds to static information of users, at a specific timestamp. For example, Ren and Tang [15] presented a vehicle location privacy protection framework based on GeoInd to compute traffic density distribution at a specific timestamp. Wang et al. [16] and Qui et al. [17] presented mobile crowdsensing frameworks in which workers' density distribution is computed in a privacypreserving manner using GeoInd. On the contrary, our work aims to estimate users' transition probabilities, which correspond to moving information of users, from perturbed trajectory datasets collected using GeoInd. To the best of our knowledge, this is the first attempt to estimate transition probability between two locations under GeoInd setting. The contributions of this paper can be briefly summarized as follows: We develop a privacy-preserving framework based on GeoInd for the computation of transition probability between two locations. In particular, the proposed scheme exploits the functionality of the expectation-maximization (EM) algorithm to precisely estimate hidden transition patterns based on perturbed trajectory dataset. Through experiments with real data, it is demonstrated that a good estimation of transition pattern can be achieved with the proposed framework.

Geo-indistinguishability
GeoInd is an extended version of differential privacy with a distance metric to provide a privacy-preserving mechanism for location data [10]. Formally, �-GeoInd is defined as follows: Definition 1 (�-GeoInd) Let X be a set of possible user locations and Y be a set of reported locations. It is commonly assumed that X is equal to Y. Let us assume that a randomized mechanism, K, probabilistically generates a perturbed location from a true location of a user. Then, K satisfies �-GeoInd, if and only if for (1) all x; x 0 2 X and (2) any output location, y 2 Y, the following equation is satisfied: where d(x, x 0 ) is the distance between x and x 0 . The parameter �, which corresponds to a privacy budget, determines a trade-off between the level of privacy and data utility. That is, smaller values of � ensure a stronger privacy guarantee, introducing larger perturbation to the true location; whereas larger values of � guarantee a weaker privacy but introduce smaller noise to the true location. This definition denotes that, given a reported location y that is an output of the randomized mechanism K, the ability of an adversary to identify whether the actual location of a user is x or x 0 is limited by the privacy budget (i.e., �) and distance between x and x 0 . This implies that the closer two locations are, the more indistinguishable they are.

Problem definition
In this subsection, we introduce key notations and formally define the problem. The main notations used in the rest of this paper are summarized in Table 1.
Let us assume that the entire area is partitioned into m × n grid G. Then, if the current location of a specific user belongs to a grid g 2 G at time t, the location of this user is represented as (g, t). The trajectory of the k-th user over r timestamps is represented as the sequence of r locations, such as TR k = {(g k1 , t k1 ), (g k2 , t k2 ), . . ., (g kr , t kr )}. Given the k-th user's trajectory TR k , let TR kp ¼ fðg 0 k1p ; t k1 Þ; ðg 0 k2p ; t k2 Þ; . . . ; ðg 0 krp ; t kr Þg be the corresponding perturbed trajectory obtained by the perturbation mechanism of GeoInd. Here, g 0 ip corresponds to the perturbed location generated from the i-th true location g i using GeoInd. Here, to differentiate between the perturbed and true locations, we use g to represent a true location and g 0 to represent a perturbed location. We will also omit the subscript, k, if it is clear from context. Let DB TR p ¼ fTR 1p ; TR 2p ; . . . TR fp g be the collection of perturbed trajectories received from all users that are maintained by the data aggregation server. Here, f is the number of users. For each TR kp ¼ fðg 0 1p ; t 1 Þ; ðg 0 2p ; t 2 Þ; . . . ; ðg 0 rp ; t r Þg in DB, we can extract a set of perturbed location pairs in two adjacent timestamps, such as TR kp pair ¼ fðg 0 Then, let DB be a collection of all such pairs of locations that can be obtained from all perturbed trajectories in DB TR p . That is, DB is formally defined as follows: where S represents the union all operator which allows duplicates. Let p(g x ! g y ) be the transition probability from g x to g y , which denotes the probability that a user currently located at g x moves to g y in the next timestamp. Given DB, a straightforward solution for obtaining transition probabilities is to directly compute p(g x ! g y ) based on the perturbed trajectories, such as where cntððg 0 x ; g 0 y Þ; DBÞ denotes the number of times that the location pair ðg 0 x ; g 0 y Þ appears in DB. However, this straightforward approach cannot accurately compute the transition probabilities because it does not consider the effect of the location perturbation mechanism of GeoInd. Thus, in this study, we aim to propose a novel scheme to accurately estimate the transition probabilities for all pairs of locations in G with the collection of perturbed trajectories.

Privacy-preserving computation of transition probability between locations
In this section, we describe our privacy-preserving framework for effectively computing transition probabilities based on the collection of perturbed trajectories collected using GeoInd.
Prob. that a user's location is g x pðg 0 y j g x Þ Prob. that g 0 y is randomly generated from g x pðg x j g 0 y Þ Prob. that a true location of g 0 y is g x https://doi.org/10.1371/journal.pone.0278744.t001

Privacy-preserving trajectory collection
There are two mechanisms to achieve �-GeoInd: the Laplace and optimization mechanisms. The Laplace mechanism is simple but it is known to introduce large noise to the collected data, resulting in a lower utility. On the contrary, the optimization mechanism can provide the maximum utility for the collected data, while satisfying �-GeoInd. In the optimization mechanism of GeoInd, the server first computes the obfuscation matrix, OM, by solving a linear programming problem. Let d(�, �) be the distance metric between two grids in G, such as the Euclidean distance or the Manhattan distance. Let us further assume that κ G be the prior probability distribution on users' possible locations. κ G can be either defined by a uniform distribution or inferred from the distribution of available historical data. Then, the obfuscation matrix OM can be obtained by solving the following linear programming problem [18,19]: min : Here, OM [u, v] represents the probability that a perturbed location g 0 v is randomly generated from a true location g u . The first constraint corresponds to the definition of GeoInd. The second constraint P g 0 v 2G OM½u; v� ¼ 1 denotes that the sum of probabilities that each location is perturbed to other locations should be 1. The third constraint OM [u, v]�0 denotes that the probability is no less than 0. It is well known that the number of constraints of the abovementioned linear optimization problem is proportional to OðjGj 3 Þ; therefore, the computational overhead of the optimal mechanism is significant. To alleviate the computational overhead of the optimal mechanism, several studies based on an approximation technique have been conducted in the literature [18,19]. After computing the obfuscation matrix OM, the server distributes it to each user. Once receiving OM, each user perturbs his/her true trajectory according to the probabilities encoded in OM. Let us focus on the k-th user's trajectory TR k = {(g 1 , t 1 ), (g 2 , t 2 ), . . ., (g r , t r )}. Given the k-th user's trajectory TR k , the corresponding perturbed trajectory TR kp ¼ fðg 0 1p ; t 1 Þ; ðg 0 2p ; t 2 Þ; . . . ; ðg 0 rp ; t r Þg is computed in a way that the perturbed location g 0 ip is randomly generated from the true location g i based on the probabilities encoded in OM. Finally, the perturbed trajectory, TR kp , is sent to the data aggregation server. Note that during this process, true locations along users' trajectories are not exposed to the outside of their devices, because the trajectory perturbation is performed on their devices.

Computation of transition probability with a collection of perturbed trajectories
In this subsection, we present the proposed method that effectively computes the transition probability between pairs of locations in G based on a collection of perturbed trajectories, which are collected under �-GeoInd. Given DB, we first compute the joint probability P(g x , g y ), which denotes the probability that a user's current location is g x and next location is g y , by leveraging the functionality of the EM algorithm. We next compute the transition probability, p(g x ! g y ), using the joint probability P(g x , g y ).
EM is an iterative algorithm that sequentially runs two steps, namely, E-step (expectation) and M-step (maximization). In the E-step, the expected value of the likelihood is computed based on the current parameters and observed variables, whereas in the M-step, an estimation on the parameters is performed to maximize the likelihood function. To leverage the EM algorithm, we first need to define a likelihood function. The collection of location pairs, DB, can be viewed as the set of observed variables, O = {o 1 , o 2 , . . ., o |DB| }, where the i-th observed variable o i corresponds to the i-th perturbed location pair ðg 0 cur i ; g 0 nex i Þ 2 DB. We also introduce the set of latent variables Z = {z 1 , z 2 , . . ., z |DB| } where z i , which is associated to o i , is the mn × mn matrix. Here, z i [u, v] is 1, if o i (i.e., the perturbed location pair ðg 0 cur i ; g 0 nex i Þ) is randomly generated from the true location (g u , g v ) by GeoInd. Otherwise, z i [u, v] is set to 0. For simplicity, we use s j,k to represent the true location pair (g j , g k ).
Given the set of observed variables O and the set of latent variables Z, the complete data likelihood function is defined as follows: Here, π j,k represents the probability that the current and next locations of a user are g j and g k respectively that we aim to estimate by using EM. That is, π j,k corresponds to P(g j , g k ). Obviously, P mn j¼1 P mn k¼1 p j;k equals to 1. Additionally, P(o i |s j,k ) corresponds to Pððg 0 cur i ; g 0 nex i Þjðg j ; g k ÞÞ. That is, P(o i |s j,k ) denotes that the perturbed location pair ðg 0 cur i ; g 0 nex i Þ is randomly generated from the true location pair (g j , g k ) by GeoInd. We note that the perturbation process of GeoInd is independently applied to each true location in a trajectory. Hence, P(o i |s j,k ) can be computed as follows: Note that by the definition of the obfuscation matrix, OM[j, cur_i] equals to Pðg 0 cur i jg j Þ. Given the Eq (5), the log-likelihood function is defined as follows: In the last line of Eq (7), P(o i |s j,k ) is substituted by Eq (6). We now define the method that computes the transition probabilities using the EM algorithm.
Initialization. In the initialized phase, the initial parameter θ (0) is defined. Determining good initial values in EM is vital to reducing the number of iterative steps until convergence and enhancing the accuracy of parameter estimation. In this paper, we use two different schemes to determine initial values: • Uniform initialization: In the first scheme, π j,k is initialized with a uniform value, such as p j;k ¼ 1 ðm�nÞ 2 , which satisfies P mn j¼1 P mn k¼1 p j;k ¼ 1.
• Distance-based initialization: Intuitively, as the distance between two locations decreases, the transition probability between them increases. Accordingly, in the second scheme, π j,k is initialized with a value that is inversely proportional to the distance between g j and g k (i.e., p j;k / 1 dðg j ;g k Þ ), while satisfying P mn j¼1 P mn k¼1 p j;k ¼ 1.

E-
Step. In this step, the conditional expectation of the latent variables Z is estimated based on the current parameter θ (h) . The E-step is stated as follows: Here, t ðhÞ i;½j;k� can be viewed as a posterior probability that given the current parameter θ (h) , the perturbed location pair ðg 0 cur i ; g 0 nex i Þ is generated from the true location pair (g j , g k Note that the last line of Eq (9) is rewritten using Eq (6).

M-
Step. This step finds the parameters θ that maximize the expectation function, Q(θ|θ (h) ), in the E-step. This task can be viewed to find the extreme value of Eq (8) with respect to the constraint P mn j¼1 P mn k¼1 p j;k ¼ 1. Therefore, we exploit the Lagrange multiplier method by defining the Lagrange function as follows: Then, the first-order partial derivative of Lðy; p; mÞ with respect to π j,k is obtained as follows: By utilizing the constraint to get the optimal value, we have the following: According to Eqs (11) and (12), the parameter that maximizes the expectation function, Q(θ|θ (h) ), is defined as follows: Thus, we can obtain the updated parameter p ðhþ1Þ j;k using Eq (13), in which t ðhÞ i;½j;k� is already computed in the previous E-step.
Computation of transition probability. The abovementioned E-step and M-step are repeated until convergence. After computing the parameters using EM, the transition probability from g x to g y is computed as follows: We note that the parameter π x,y is already estimated in the previous E-Step and M-Step.

Experiments
In this section, we present the experiments we carried out to evaluate the proposed approach. First we describe the experimental setup and then we will discuss the experimental results.

Experimental setup
We evaluate the proposed method using the T-Drive dataset [20] that contains one-week trajectories of Beijing taxis. We first extract 85707 trajectories, with a length of 10, from the T-Drive dataset. Then, the entire geographic region, represented by longitude and latitude information, is segmented into three different grids: 10 × 10, 15 × 15, and 20 × 20. Each location in trajectories is assigned to a segmented region (i.e., grid) to which it geographically belongs. In the experiments, results are reported for the following alternatives, the straightforward approach (SA), the proposed EM-based approach with the uniform initialization (EM uni ), and the proposed EM-based approach with the distance-based initialization (EM dis ). Furthermore, we compare our proposed method with the particle filter-based approach (PF) which has been extensively used in trajectory detection [21,22]. For our evaluations, we adapted the underlying particle filter-based approach technique for trajectory detection to infer users' true locations from perturbed locations. Especially, at each iteration of the particle filter, the weights of particles are updated by using the probabilities embedded in the obfuscation matrix OM of GeoInd.
To compare these four schemes, we use the mean absolute error (MAE), which is defined as follows: Here, p est (g x ! g y ) is the transition probability estimated based on the perturbed trajectory dataset, whereas p(g x ! g y ) is the true transition probability computed based on the actual trajectory dataset.

Fig 2 illustrates MAE values versus a varying privacy budget �.
In the experiments, � varies from 0.5 to 2.0, while the size of grid is fixed to 10 × 10. Key observations based on Fig 2 can be summarized as follows: As expected, the error rate increases as the privacy budget decreases. This is because a smaller � value provides stronger privacy. However, in terms of utility, a smaller � results in a lower utility of the collected trajectory dataset. This in turn leads to a decreased estimation accuracy when computing the transition probability. On the contrary, a larger value of � provides a weaker privacy guarantee, while introducing less perturbation to true locations. This in turn leads to an increased estimation accuracy of the transition probability.
Among the four alternatives, the proposed EM-based schemes, EM uni and EM dis , exhibit the better performance compared to the straightforward approach SA and the particle filterbased approach PF. Especially, the performance gain of the proposed EM-based schemes over SA and PF is more pronounced at a smaller � value, which provides stronger privacy. Considering that many real-world applications require stronger privacy protection guarantees against attackers with arbitrary backgrounds, these results indicate that the proposed method is more practical for real-world applications. Between the two EM-based schemes, EM dis slightly outperforms EM uni at all privacy bugets, which indicates that the distance-based initialization scheme is more promising than a uniform value scheme.
To further investigate the validity of the experimental results, in Fig 3, we plot a heatmap for the estimated transition probability distribution of each method with � being set to 1.0. The size of grid is set to 10 × 10. For comparison purposes, we also plot the true transition probability distribution that is obtained with the actual trajectory dataset. As observed in the figure, the  probability distribution obtained by SA and PF is quite dissimilar to the true distribution. On the contrary, with the proposed EM-based schemes, we can obtain the estimated probability distributions that are highly similar to the true one. Between the two EM-based schemes, the estimated probability distributions obtained by EM dis are more similar with the true distributions than those computed by EM uni . The experimental results presented in Figs 2 and 3 indicate that the proposed EM-based approach can achieve higher precision in the computation of the transition probability with trajectory dataset collected in a privacy-preserving manner using GeoInd than that of the straightforward and particle filter-based approaches. Fig 4 shows MAE values with respect to various grid sizes. In the experiment, three different grid sizes are used with � fixed to 0.5. Similar experimental results can be observed with various grid sizes. As shown in Fig 4, for all grid sizes, the proposed EM uni and EM dis significantly outperform SA and PF regarding MAE, verifying the robustness of the proposed method against the grid size. Furthermore, the EM scheme using the distance-based initialization achieves slightly better performance than that using the uniform-based initialization. Finally, we validate the convergence of the proposed EM-based schemes. In Fig 5, we plot the MAE values versus the number of iterations in the EM process. In the experiments, � is fixed to 1.0, while three different grid sizes are used. As shown in the figure, we can observe a significant drop in the MAE values during the first a few iteration. Beyond this point, the MAE values stabilize. These experiment results verify that for the proposed EM-based scheme that computes the transition probability, the number of iterations guaranteeing convergence is small. This in turn indicates that the additional computational overhead of the proposed scheme incurred when using the EM algorithm is not significant.

Related work
GeoInd have been used in diverse application domains. Here, we present some application areas of GeoInd. Vehicle networks have recently emerged as a promising solution to improve driving experiences and road safety. In [23], Zhou et al. proposed edge-assisted vehicle networks for improving the service quality in the vehicle networks where GeoInd is deployed at the edge nodes to protect the true location of the vehicle. In the proposed system, a vehicle first submits a service request along with its location to an edge node at which the GeoInd-based mechanism is executed to protect the actual location of the vehicle, and the service request with the perturbed vehicle's location is then forwarded to the service provider which returns the required service to the requesting vehicle. To protect location privacy in location-aware social networks, GeoInd combined with homomorphic encryption is used for privacy-preserving nearby friend discovery [24]. Spatial crowdsourcing is a platform where individual users are engaged to collect, analyze, and disseminate their surrounding information. Wang et al. [16] proposed a differential geo-obfuscation mechanism to protect the workers' true location during task assignment by an spatial crowdsourcing platform. Yan et al. [25] developed a spatial crowdsourcing framework which can protect the privacy of the workers' trajectories. In the developed framework, the GeoInd mechanism is used to protect a worker's shortest path from the source to destination. Qui et al. [17] investigated location-privacy protection in a vehicle-based spatial crowdsourcing framework using GeoInd. To address the location privacy issues raised in ride-sharing applications such as Uber, Waze, and Lyft, Tong et al. [26] proposed a scheduling scheme which exploits GeoInd to protect the location information of ridesharing users. We note that the existing works fucus on either protecting users' location privacy using GeoInd or extracting static information, such as population density distribution, from perturbed location datasets collected under GeoInd setting. On the contrary, the objective of this work is to estimate users' moving information, such as transition patterns, with perturbed trajectory datasets collected using GeoInd.
There have been extensive studies to leverage the concept of differential privacy for publishing trajectory data in a privacy-preserving manner. Hua et al. [27] proposed a differential-privacy-based scheme for publishing time-serial trajectory data. The proposed scheme leverages an exponential mechanism to probabilistically cluster locations based on their distance, and then relies on the Laplace mechanism to add a random noise to the count of trajectories in a cluster. DP-Star [28] is a differential-privacy-based framework for publishing trajectory data with strong utility. DP-Star generates synthetic trajectories that satisfy �-differential privacy, while maintaining high utility. SafePath [29], which is a privacy-preserving algorithm for publishing trajectories, structures trajectories as a noisy prefix tree and publishes differentially-private trajectories, while retaining data utility. Ou et al. [30] proposed two lagrange multiplier- based differentially private approaches, UD-LMDP and UC-LMDP, to address privacy issues arising when publishing mutually correlated trajectories. The authors introduced an n-body Laplace framework which aims to prevent adversaries from inferring a social relation from the mutual correlation between two users' trajectories. Chen et al. [11] proposed RNN-DP, which is a differential privacy scheme based on a recurrent neural network, to protect the privacy of real-time trajectory data. RNN-DP exploits a recurrent neural network to efficiently predict trajectory, while protecting the location privacy of users using differential privacy. OPTDP [31] is an optimal personalized trajectory differential privacy mechanism for trajectory data publishing. Zhao et al. [32] introduced a method to build the SR tree, which is based on R-tree, using the trajectory sequence, and then adds random noise to the nodes in the SR tree. Then, the noise SR-Tree indexes are used to answer user's spatial queries. Liu et al. [33] developed a trajectory data publication scheme that exploits the staircase mechanism of differential privacy.
Differential privacy has been also used to process and analyze trajectory data in diverse application areas. PLDP-TD [34] supports personalized-location differentially private data analysis on trajectory databases. To answer queries in a differentially private way, PLDP-TD builds a personalized noisy trajectory tree which stores sub-trajectories of a trajectory database along with their privacy protection levels. That is, a personalized privacy level is assigned to each node of the tree depending on the privacy protection requirements of locations. Kim et al. [35] leverages local differential privacy, which is a localized version of differential privacy, to collect indoor positioning data from users, while protecting their location privacy. Con-Crowd-DP [36] is a differentially private framework for mobile crowdsourcing applications in which the mobile users upload the location-related task results to the server for obtaining the rewards. In ConCrowd-DP, to protect of participating users' location privacy, perturbed locations are used instead of users' true locations, when reporting the location-related task results to the server. Deldar and Abadi [37] proposed DPLG to efficiently and accurately answer spatial queries on moving objects databases. DPLG is based on the combination of differential privacy and location generalization which enable to efficiently and accurately process spatial queries by reducing the number of locations and minimizing query errors.

Conclusion
The most common analysis task in applications that utilize mobility information of users is to compute the transition probability between two consecutive timestamps. Computing transition probabilities from a collection of privacy-preserved trajectories is challenging, because of a loss in the utility of collected dataset caused by a privacy-preserving mechanism. To address this challenge, in this paper, we proposed a novel scheme to compute the transition probability with the collection of perturbed trajectories collected using GeoInd. The proposed method leveraged the EM algorithm to precisely estimate hidden transition patterns based on a perturbed trajectory dataset. Experimental results with real datasets verified that a good estimation of transition pattern can be achieved with the proposed method.